Data description

Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.

me5 <- load_object("model-2019-04-09-19-02-16.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me5)
## Summary for model '/var/folders/2z/0sxx_8ts2pxcp028sy41wfgc0000gn/T//RtmpeNqPAL/file1219c690e12c'
## Saved parameters: B Rho EnvRho Tau 
## MCMC ran in parallel for 34.541 minutes at time 2019-04-09 18:27:42.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#me5$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(me5$R$B),4)))
## Maximum Rhat value for Beta: 1.136
cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(me5$R$Rho),4)))
## Maximum Rhat value for Rho: NA
cat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(me5$R$EnvRho),4)))
## Maximum Rhat value for EnvRho: 1.1004
cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(me5$R$Tau),4)))
## Maximum Rhat value for Tau: 1.4694
me5$n.eff
## $B
##      [,1] [,2] [,3]
## [1,]   29  155   32
## [2,]  257  547  283
## [3,]  756 1000 1000
## [4,]  103  593   94
## [5,]  181  204  177
## 
## $Rho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1   86   66  284  232
## [2,]   86    1  132  204  350
## [3,]   66  132    1  195   53
## [4,]  284  204  195    1  195
## [5,]  232  350   53  195    1
## 
## $EnvRho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1   54   39   53   48
## [2,]   54    1  231  129  277
## [3,]   39  231    1  189  273
## [4,]   53  129  189    1  181
## [5,]   48  277  273  181    1
## 
## $Tau
##      [,1] [,2] [,3] [,4] [,5]
## [1,]   12  167  129  116  376
## [2,]  167   29  411   82  216
## [3,]  129  411   33  491  159
## [4,]  116   82  491   21  251
## [5,]  376  216  159  251   59
me5$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Environment filtering 5 species

JSDM

Model for \(\textbf{EnvEvenSp5}\)

data<-sim_data$EnvEvenSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

Y_cor<-cor(data$Y)

to_prec<-function(m){
  n<-dim(m)[1]
  Tau_n<-matrix(nrow=n, ncol=n)
  for (j in 1:n) {
   for (k in 1:n){
    Tau_n[j, k] <-  -m[j, k]/sqrt((m[j,j]*m[k,k]))
   }
  }
  return(Tau_n)
}

#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
title("Tau signif")

GJAM

##to setup chains parameters

data<-sim_data$EnvEvenSp5
#fit_gjam(data,2000,1000,"./gjam_models/gjam5env.rda")
load_gjam(data,name="./gjam_models/gjam5env.rda")
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       9.5 3.13   5.08   16.4
## env2     31.9 7.93  20.50   50.6
## 
## Coefficient matrix B:
##              sp01   sp02    sp03   sp04   sp05
## intercept -0.0312  0.646  1.2500  0.704 -1.200
## env       -2.4200 -1.800  0.0438  1.710  2.870
## env2       0.9440 -0.527 -0.9410 -0.598 -0.665
## RMSPE      0.2780  0.327  0.3760  0.309  0.327
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025 CI_975 sig95
## sp01_intercept  -0.0312 0.0775 -0.1940  0.107      
## sp01_env        -2.4200 0.1170 -2.6500 -2.200     *
## sp01_env2        0.9440 0.0938  0.7780  1.130     *
## sp02_intercept   0.6460 0.0651  0.5220  0.771     *
## sp02_env        -1.8000 0.0911 -1.9900 -1.630     *
## sp02_env2       -0.5270 0.0741 -0.6700 -0.389     *
## sp03_intercept   1.2500 0.0672  1.1100  1.380     *
## sp03_env         0.0438 0.0465 -0.0447  0.132      
## sp03_env2       -0.9410 0.0764 -1.0900 -0.781     *
## sp04_intercept   0.7040 0.0865  0.5590  0.881     *
## sp04_env         1.7100 0.0897  1.5500  1.890     *
## sp04_env2       -0.5980 0.0608 -0.7170 -0.486     *
## sp05_intercept  -1.2000 0.0953 -1.3700 -1.000     *
## sp05_env         2.8700 0.1380  2.6300  3.170     *
## sp05_env2       -0.6650 0.0610 -0.7840 -0.545     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.325, and the DIC is 25294.  Computation involved 2000 Gibbs steps, with a burnin of 1000.

HMSC

#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="F_t",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
  if (label=="F_t"){
    Y_data = subset(data, select = -env)
    ns<- ncol(Y_data)
    np <- nrow(Y_data)
    X<-scale(poly(data$env[1:np], 2))
    colnames(X)<-c("env","env2")
    studyDesign = data.frame(sample = as.factor(1:np))
    rL = HmscRandomLevel(units = studyDesign$sample)
    m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
             studyDesign = studyDesign, ranLevels = list(sample = rL))
    m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
    save(m, file = name)
    return(m)
  }
  if (label=="L_d"){
    return(load_object(name))
  }
}

data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")

Convergence:

hm_conv<-function(mod){
  codaList = convertToCodaObject(mod)

  #convergence histograms
  hist(effectiveSize(codaList$Beta), main="ess(beta)")
  hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf, main="psrf(beta)")
  
  hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)")
  hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf, main="psrf(omega)")
}

hm_conv(hm_mod)

Study of interactions

hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
  getOmega = function(a,r=1)
  return(crossprod(a$Lambda[[r]]))
  ns<-mod$ns
  postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
  postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
  
  postOmega<-abind(postOmega1,postOmega2,along=3)
  postOmegaMean = apply(postOmega,c(1,2),mean)
  postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
  postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
  
  postR<-array(dim=c(ns,ns,nchains*nsamples))
  for(i in 1:dim(postOmega)[3])
  postR[,,i]<-stats::cov2cor(postOmega[,,i])
  
  
  postRMean = apply(postR,c(1,2),mean)
  postRUp=apply(postR,c(1,2),quantile,0.95)
  postRLo=apply(postR,c(1,2),quantile,0.05)
  
  Tau = solve(postOmegaMean)
  Tau_n = cov2cor(Tau)
  
  
  Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
  
  # Omegacor<- computeAssociations(m)
  # supportLevel<- 0.95
  # toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
  # corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
  
  par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  corrplot(cor(hm_mod$Y), diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Correlation cor(Y)")
  corrplot(postRMean, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("R")
  corrplot(Toplot_R, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Plot only non zero value")
  corrplot(Tau_n, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Partial correlation matrix")
  corrplot(interact, diag = FALSE, order = "FPC",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("True interactions")
}


hm_inter(hm_mod)

Environmental filtering 10 species

JSDM

me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)
## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#me10$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(me10$R$B),4)))
## Maximum Rhat value for Beta: 1.5602
cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(me10$R$Rho),4)))
## Maximum Rhat value for Rho: NA
cat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(me10$R$EnvRho),4)))
## Maximum Rhat value for EnvRho: 1.7953
cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(me10$R$Tau),4)))
## Warning in max(me10$R$Tau): no non-missing arguments to max; returning -Inf
## Maximum Rhat value for Tau: -Inf
me10$n.eff
## $B
##       [,1] [,2] [,3]
##  [1,]   34   51   38
##  [2,]   33  121   31
##  [3,]  130  137   82
##  [4,]  271 1000  645
##  [5,]  824  593 1000
##  [6,]  756  670 1000
##  [7,]  729  670  259
##  [8,]   58   69   61
##  [9,]   45   47   47
## [10,]   11   47   10
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1  119   49   54   84  102   66  519  121   186
##  [2,]  119    1  162  121   81  136  184  283  224   219
##  [3,]   49  162    1  755   53   34  120   72  256   217
##  [4,]   54  121  755    1  357   80   82   75   22   273
##  [5,]   84   81   53  357    1  269   76   84   62    59
##  [6,]  102  136   34   80  269    1  249  404  136    78
##  [7,]   66  184  120   82   76  249    1  391  363    30
##  [8,]  519  283   72   75   84  404  391    1   36    66
##  [9,]  121  224  256   22   62  136  363   36    1   157
## [10,]  186  219  217  273   59   78   30   66  157     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1   60   74   80   59   39   39   34   89    18
##  [2,]   60    1   36   56   44   49   39   86  230     9
##  [3,]   74   36    1  144  134  112  110   51  249    23
##  [4,]   80   56  144    1  505  379  214   67   74    14
##  [5,]   59   44  134  505    1  632 1000  113  106    16
##  [6,]   39   49  112  379  632    1  801  103   86    14
##  [7,]   39   39  110  214 1000  801    1  128  102    14
##  [8,]   34   86   51   67  113  103  128    1   42    17
##  [9,]   89  230  249   74  106   86  102   42    1    23
## [10,]   18    9   23   14   16   14   14   17   23     1
me10$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
data<-sim_data$EnvEvenSp10
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
  par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
  corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Correlation cor(Y)")
  corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("EnvRho")
  corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("EnvRho signif")
  corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Rho")
  corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Rho signif")
  corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("True interactions")
  #corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  #title("Tau")
  #corrplot(Tau_n*(!me10$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  #title("Tau")
}

plot_cor_jsdm(me10,data$Y)

Gjam

data<-sim_data$EnvEvenSp10

#fit_gjam(data,5000,500,"./gjam_models/gjam10env.rda")
load_gjam(data,name="./gjam_models/gjam10env.rda")
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env     233.0 58.0    120  354.0
## env2     56.7 11.8     37   82.9
## 
## Coefficient matrix B:
##              sp01    sp02   sp03   sp04   sp05   sp06   sp07   sp08  sp09
## intercept -0.2060  0.3800  1.310  0.889  1.020  0.742  0.638  0.140 0.332
## env       -2.0500 -2.5700 -2.300 -1.030 -0.189  0.542  1.090  1.740 2.520
## env2       0.0656 -0.0531  0.338 -0.533 -0.805 -0.838 -0.690 -0.425 0.640
## RMSPE      0.3300  0.2980  0.317  0.371  0.414  0.421  0.374  0.345 0.298
##             sp10
## intercept -0.577
## env        2.370
## env2       0.677
## RMSPE      0.276
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept  -0.2060 0.0670 -0.33700 -0.0733     *
## sp01_env        -2.0500 0.1190 -2.29000 -1.8300     *
## sp01_env2        0.0656 0.0992 -0.14100  0.2380      
## sp02_intercept   0.3800 0.1030  0.18000  0.5510     *
## sp02_env        -2.5700 0.1370 -2.83000 -2.2900     *
## sp02_env2       -0.0531 0.0790 -0.19700  0.1120      
## sp03_intercept   1.3100 0.0829  1.15000  1.4700     *
## sp03_env        -2.3000 0.1260 -2.55000 -2.0600     *
## sp03_env2        0.3380 0.0712  0.20500  0.4860     *
## sp04_intercept   0.8890 0.0585  0.77400  1.0000     *
## sp04_env        -1.0300 0.0659 -1.16000 -0.9070     *
## sp04_env2       -0.5330 0.0580 -0.64200 -0.4160     *
## sp05_intercept   1.0200 0.0621  0.90800  1.1500     *
## sp05_env        -0.1890 0.0488 -0.28500 -0.0949     *
## sp05_env2       -0.8050 0.0759 -0.96000 -0.6590     *
## sp06_intercept   0.7420 0.0635  0.61600  0.8650     *
## sp06_env         0.5420 0.0558  0.43200  0.6460     *
## sp06_env2       -0.8380 0.0634 -0.96100 -0.7130     *
## sp07_intercept   0.6380 0.0795  0.50600  0.8160     *
## sp07_env         1.0900 0.0805  0.93900  1.2500     *
## sp07_env2       -0.6900 0.0595 -0.80300 -0.5700     *
## sp08_intercept   0.1400 0.0759  0.00798  0.2950     *
## sp08_env         1.7400 0.1240  1.48000  1.9600     *
## sp08_env2       -0.4250 0.1000 -0.65400 -0.2750     *
## sp09_intercept   0.3320 0.0539  0.22600  0.4370     *
## sp09_env         2.5200 0.1660  2.22000  2.8500     *
## sp09_env2        0.6400 0.0700  0.50400  0.7780     *
## sp10_intercept  -0.5770 0.0696 -0.70200 -0.4380     *
## sp10_env         2.3700 0.1130  2.16000  2.6000     *
## sp10_env2        0.6770 0.0820  0.49900  0.8200     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.348, and the DIC is 114934.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gje10<-load_object("./gjam_models/gjam10env.rda")

#to check posterior density of s in Sigma 
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))

HMSC

data<-sim_data$EnvEvenSp10
#hm_mod<-fit_hmsc(data,"F_t",name="./HMmodels/hm10env.rda" )
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm10env.rda" )

hm_conv(hm_mod)

hm_inter(hm_mod)

Environmental filtering 20 species

JSDM

me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)
## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#me20$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(me20$R$B),4)))
## Maximum Rhat value for Beta: 1.2174
cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(me20$R$Rho),4)))
## Maximum Rhat value for Rho: NA
cat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(me20$R$EnvRho),4)))
## Maximum Rhat value for EnvRho: 1.3236
cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(me20$R$Tau),4)))
## Warning in max(me20$R$Tau): no non-missing arguments to max; returning -Inf
## Maximum Rhat value for Tau: -Inf
me20$n.eff
## $B
##       [,1] [,2] [,3]
##  [1,]  114  142  148
##  [2,]   39   89   30
##  [3,]   21   99   21
##  [4,]   23   30   22
##  [5,]   26   46   30
##  [6,]   88  141   85
##  [7,]   39   35   42
##  [8,]  675  220   85
##  [9,] 1000  265 1000
## [10,]  674  700  218
## [11,]  598  276 1000
## [12,]  180  413  334
## [13,]  438  378  386
## [14,]  366  129  155
## [15,]   56   44   83
## [16,]   53   97   63
## [17,]   23   45   26
## [18,]   42   59   30
## [19,]   19   68   21
## [20,]   37  116   30
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1  181  118  141   73   46   71   89  138    25    49    63   247
##  [2,]  181    1   78  200  111   18  147  186   58    40   121    74    37
##  [3,]  118   78    1   76   43   28   20   59  124   212    85    91    74
##  [4,]  141  200   76    1  354  143  269  122  173    42    66    90    29
##  [5,]   73  111   43  354    1  533  128  100  189   370   176    57    42
##  [6,]   46   18   28  143  533    1  661   73  175    97   963    78    68
##  [7,]   71  147   20  269  128  661    1  212  143   324   100    84   125
##  [8,]   89  186   59  122  100   73  212    1  107   395   156   115   115
##  [9,]  138   58  124  173  189  175  143  107    1   813   131   143   109
## [10,]   25   40  212   42  370   97  324  395  813     1   106    67    49
## [11,]   49  121   85   66  176  963  100  156  131   106     1   431   244
## [12,]   63   74   91   90   57   78   84  115  143    67   431     1   145
## [13,]  247   37   74   29   42   68  125  115  109    49   244   145     1
## [14,]  567  119   24  413   41  348   26  353  197   135   121   178   318
## [15,]  111  174   59  353   54   31   26  295   42    48  1000    66    75
## [16,]   77  136   59   56   76   47  402   80   44   101   245   450    53
## [17,]  105   82   67   58   70  171   44   80   51    43    57    88    30
## [18,]  316  422   58   81  425  152  110  131   55   164   143    72    30
## [19,]   63  212  130  112  125  110  199   78   56    29    81    45    31
## [20,]  143  364  127  489  821  208  563   80  273    74   318   679    58
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]   567   111    77   105   316    63   143
##  [2,]   119   174   136    82   422   212   364
##  [3,]    24    59    59    67    58   130   127
##  [4,]   413   353    56    58    81   112   489
##  [5,]    41    54    76    70   425   125   821
##  [6,]   348    31    47   171   152   110   208
##  [7,]    26    26   402    44   110   199   563
##  [8,]   353   295    80    80   131    78    80
##  [9,]   197    42    44    51    55    56   273
## [10,]   135    48   101    43   164    29    74
## [11,]   121  1000   245    57   143    81   318
## [12,]   178    66   450    88    72    45   679
## [13,]   318    75    53    30    30    31    58
## [14,]     1   111   183    58    75    52    86
## [15,]   111     1   221   274    60   333    59
## [16,]   183   221     1   128    70   469   101
## [17,]    58   274   128     1   233   118   159
## [18,]    75    60    70   233     1   117   161
## [19,]    52   333   469   118   117     1   172
## [20,]    86    59   101   159   161   172     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1   99   33   72   38  483  109   72  139   273   137   356   205
##  [2,]   99    1   68  111   59   42   33   28   34    36    43    44    50
##  [3,]   33   68    1   51   48   33   47   31   30    32    41    28    39
##  [4,]   72  111   51    1   36   51   22   27   40    40    60    49    58
##  [5,]   38   59   48   36    1   32   55   57   50    35    66    42    69
##  [6,]  483   42   33   51   32    1   58   70  102   233    76   227   133
##  [7,]  109   33   47   22   55   58    1   65   61    50    67    62    74
##  [8,]   72   28   31   27   57   70   65    1  228   224   195   178   357
##  [9,]  139   34   30   40   50  102   61  228    1   351   520   689  1000
## [10,]  273   36   32   40   35  233   50  224  351     1   169   794   489
## [11,]  137   43   41   60   66   76   67  195  520   169     1   310  1000
## [12,]  356   44   28   49   42  227   62  178  689   794   310     1   856
## [13,]  205   50   39   58   69  133   74  357 1000   489  1000   856     1
## [14,]  312   74   32   37   48  149   91  167  351   214   287   456   201
## [15,]   62   38   43   94   58   90   68  152  184   131   103   114   112
## [16,]  148   62   27  315   28  174   39   56  128   121    76   211   132
## [17,]   58   59   56   30   32   90   20   29   36    56    29    47    48
## [18,]   44   34   25   53   64   42   34   52   77    39   123    57   108
## [19,]  100   37  127   42   53   33   50   23   25    29    30    32    31
## [20,]   46   18   15   29   20   89   35   34   56    47    45    73    68
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]   312    62   148    58    44   100    46
##  [2,]    74    38    62    59    34    37    18
##  [3,]    32    43    27    56    25   127    15
##  [4,]    37    94   315    30    53    42    29
##  [5,]    48    58    28    32    64    53    20
##  [6,]   149    90   174    90    42    33    89
##  [7,]    91    68    39    20    34    50    35
##  [8,]   167   152    56    29    52    23    34
##  [9,]   351   184   128    36    77    25    56
## [10,]   214   131   121    56    39    29    47
## [11,]   287   103    76    29   123    30    45
## [12,]   456   114   211    47    57    32    73
## [13,]   201   112   132    48   108    31    68
## [14,]     1    80   101    40    47    33    48
## [15,]    80     1  1000    63   127    19   172
## [16,]   101  1000     1    96   337    27   188
## [17,]    40    63    96     1    33   355    35
## [18,]    47   127   337    33     1    17    42
## [19,]    33    19    27   355    17     1    18
## [20,]    48   172   188    35    42    18     1
me20$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
data<-sim_data$EnvEvenSp20
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))


plot_cor_jsdm(me20,data$Y)

Gjam

data<-sim_data$EnvEvenSp20

#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

load_gjam(data,name="./gjam_models/gjam20env.rda")
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       335 61.7    222    461
## env2      157 28.3    110    219
## 
## Coefficient matrix B:
##             sp01    sp02   sp03   sp04   sp05   sp06   sp07   sp08   sp09
## intercept -0.533  0.0688  0.500  0.144  0.327  0.498  0.720  0.896  0.906
## env       -1.850 -2.5200 -2.480 -2.020 -1.820 -1.760 -1.220 -0.753 -0.516
## env2       0.187  0.3400  0.451 -0.279 -0.507 -0.493 -0.651 -0.641 -0.815
## RMSPE      0.337  0.3000  0.309  0.331  0.332  0.333  0.365  0.398  0.406
##             sp10   sp11   sp12   sp13   sp14   sp15   sp16   sp17     sp18
## intercept  1.130  1.120  1.050  0.689  0.610  0.435  0.440  0.113 -0.07100
## env       -0.200  0.133  0.551  0.870  1.140  1.730  2.460  1.990  2.35000
## env2      -0.931 -0.855 -0.860 -0.794 -0.671 -0.452 -0.227 -0.271  0.00271
## RMSPE      0.398  0.411  0.386  0.399  0.363  0.339  0.308  0.330  0.30300
##              sp19   sp20
## intercept -0.0697 -0.329
## env        2.5800  2.190
## env2       0.3120  0.530
## RMSPE      0.3030  0.318
## 
## Coefficient matrix B:
##                Estimate     SE    CI_025  CI_975 sig95
## sp01_intercept -0.53300 0.1340 -0.764000 -0.3250     *
## sp01_env       -1.85000 0.1340 -2.100000 -1.5800     *
## sp01_env2       0.18700 0.0717  0.055600  0.3320     *
## sp02_intercept  0.06880 0.0972 -0.095900  0.2760      
## sp02_env       -2.52000 0.1560 -2.810000 -2.2200     *
## sp02_env2       0.34000 0.0720  0.206000  0.4840     *
## sp03_intercept  0.50000 0.1040  0.295000  0.6900     *
## sp03_env       -2.48000 0.1360 -2.770000 -2.2400     *
## sp03_env2       0.45100 0.0813  0.310000  0.6170     *
## sp04_intercept  0.14400 0.0783  0.009940  0.3210     *
## sp04_env       -2.02000 0.0913 -2.200000 -1.8400     *
## sp04_env2      -0.27900 0.0825 -0.432000 -0.1160     *
## sp05_intercept  0.32700 0.0843  0.172000  0.4850     *
## sp05_env       -1.82000 0.1010 -2.020000 -1.6200     *
## sp05_env2      -0.50700 0.0716 -0.653000 -0.3720     *
## sp06_intercept  0.49800 0.0786  0.342000  0.6400     *
## sp06_env       -1.76000 0.0786 -1.920000 -1.6100     *
## sp06_env2      -0.49300 0.0669 -0.631000 -0.3680     *
## sp07_intercept  0.72000 0.0628  0.598000  0.8430     *
## sp07_env       -1.22000 0.0844 -1.370000 -1.0500     *
## sp07_env2      -0.65100 0.0816 -0.819000 -0.4980     *
## sp08_intercept  0.89600 0.0702  0.766000  1.0400     *
## sp08_env       -0.75300 0.0678 -0.886000 -0.6230     *
## sp08_env2      -0.64100 0.0614 -0.759000 -0.5180     *
## sp09_intercept  0.90600 0.0712  0.776000  1.0500     *
## sp09_env       -0.51600 0.0835 -0.663000 -0.3500     *
## sp09_env2      -0.81500 0.0669 -0.942000 -0.6790     *
## sp10_intercept  1.13000 0.0655  1.000000  1.2600     *
## sp10_env       -0.20000 0.0573 -0.305000 -0.0862     *
## sp10_env2      -0.93100 0.0676 -1.060000 -0.8010     *
## sp11_intercept  1.12000 0.0684  0.987000  1.2500     *
## sp11_env        0.13300 0.0565  0.022800  0.2410     *
## sp11_env2      -0.85500 0.0630 -0.978000 -0.7310     *
## sp12_intercept  1.05000 0.0730  0.906000  1.2000     *
## sp12_env        0.55100 0.0747  0.420000  0.7010     *
## sp12_env2      -0.86000 0.0752 -1.010000 -0.7180     *
## sp13_intercept  0.68900 0.0591  0.570000  0.8000     *
## sp13_env        0.87000 0.0813  0.709000  1.0300     *
## sp13_env2      -0.79400 0.0668 -0.922000 -0.6560     *
## sp14_intercept  0.61000 0.0641  0.481000  0.7300     *
## sp14_env        1.14000 0.0675  1.020000  1.2800     *
## sp14_env2      -0.67100 0.0665 -0.800000 -0.5470     *
## sp15_intercept  0.43500 0.0588  0.325000  0.5530     *
## sp15_env        1.73000 0.1250  1.520000  1.9800     *
## sp15_env2      -0.45200 0.0641 -0.579000 -0.3340     *
## sp16_intercept  0.44000 0.0521  0.337000  0.5420     *
## sp16_env        2.46000 0.1050  2.260000  2.6700     *
## sp16_env2      -0.22700 0.0501 -0.323000 -0.1280     *
## sp17_intercept  0.11300 0.0563  0.000618  0.2230     *
## sp17_env        1.99000 0.1310  1.760000  2.2600     *
## sp17_env2      -0.27100 0.0645 -0.394000 -0.1480     *
## sp18_intercept -0.07100 0.0549 -0.181000  0.0357      
## sp18_env        2.35000 0.1040  2.140000  2.5500     *
## sp18_env2       0.00271 0.0659 -0.130000  0.1260      
## sp19_intercept -0.06970 0.0875 -0.241000  0.0747      
## sp19_env        2.58000 0.1220  2.350000  2.8300     *
## sp19_env2       0.31200 0.0517  0.216000  0.4120     *
## sp20_intercept -0.32900 0.0556 -0.433000 -0.2160     *
## sp20_env        2.19000 0.0953  2.010000  2.3800     *
## sp20_env2       0.53000 0.0616  0.409000  0.6490     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.351, and the DIC is 371075.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gje20<-load_object("./gjam_models/gjam20env.rda")

#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))

Gjam dimension reduction

data<-sim_data$EnvEvenSp20



#fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")

#load_gjam(data,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")

#to check posterior density of s in Sigma 
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
 data <- list(
    Y = subset(data, select = -env),
    X = cbind(1, scale(poly(data$env, 2))),
    covx = cov(cbind(1, scale(poly(data$env, 2)))),
    K = 3,
    J = ncol(data) - 1,
    n = nrow(data),
    I = diag(ncol(data) - 1),
    df = ncol(data)
  )
  xdata<-as.data.frame(data$X[,-1])
  colnames(xdata)<- c("env","env2")
  ydata<-as.data.frame(data$Y)
  ###formula
  rl   <- list(r = 8, N = 20)
  formula<-as.formula( ~env+ env2)
  ml  <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
  ####fit
  
  
  mod_gjam1  <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
## 
## Observations and responses:
## [1] 500  20
## Warning in .setupReduct(modelList, S, Q, n): dimension reduction requires
## reductList$N < no. responses
## 
## Dimension reduced from 20 X 20 -> 20 X 8 responses
## ===========================================================================
## expanding covariance chains
## ===========================================================================
  save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
  summary(mod_gjam1)
## 
## Sensitivity by predictor variables f:
##      Estimate     SE CI_025 CI_975
## env     0.963 0.0606  0.844   1.08
## env2    1.310 0.0698  1.180   1.45
## 
## Coefficient matrix B:
##             sp01    sp02     sp03   sp04   sp05   sp06   sp07   sp08
## intercept -0.344 -0.1010  0.00677  0.111  0.292  0.342  0.469  0.483
## env       -0.473 -0.4810 -0.48100 -0.475 -0.470 -0.463 -0.430 -0.356
## env2       0.276  0.0944  0.00263 -0.110 -0.262 -0.329 -0.441 -0.460
## RMSPE      0.406  0.4050  0.40800  0.410  0.403  0.402  0.408  0.424
##             sp09   sp10   sp11   sp12   sp13   sp14   sp15   sp16   sp17
## intercept  0.484  0.488  0.488  0.487  0.478  0.475  0.392  0.276  0.177
## env       -0.362 -0.091  0.198  0.462  0.423  0.457  0.466  0.476  0.478
## env2      -0.478 -0.484 -0.475 -0.482 -0.478 -0.454 -0.347 -0.259 -0.171
## RMSPE      0.431  0.441  0.446  0.427  0.425  0.405  0.406  0.407  0.408
##                sp18    sp19   sp20
## intercept -0.000261 -0.0922 -0.277
## env        0.481000  0.4800  0.471
## env2       0.009230  0.1010  0.246
## RMSPE      0.404000  0.4060  0.403
## 
## Coefficient matrix B:
##                 Estimate     SE  CI_025  CI_975 sig95
## sp01_intercept -0.344000 0.0705 -0.4720 -0.2020     *
## sp01_env       -0.473000 0.0296 -0.4990 -0.3980     *
## sp01_env2       0.276000 0.0687  0.1430  0.4020     *
## sp02_intercept -0.101000 0.0732 -0.2410  0.0477      
## sp02_env       -0.481000 0.0205 -0.5000 -0.4260     *
## sp02_env2       0.094400 0.0682 -0.0350  0.2330      
## sp03_intercept  0.006770 0.0748 -0.1410  0.1520      
## sp03_env       -0.481000 0.0199 -0.5000 -0.4280     *
## sp03_env2       0.002630 0.0679 -0.1350  0.1300      
## sp04_intercept  0.111000 0.0713 -0.0349  0.2420      
## sp04_env       -0.475000 0.0245 -0.4990 -0.4100     *
## sp04_env2      -0.110000 0.0675 -0.2400  0.0263      
## sp05_intercept  0.292000 0.0705  0.1490  0.4300     *
## sp05_env       -0.470000 0.0291 -0.4990 -0.3920     *
## sp05_env2      -0.262000 0.0691 -0.3960 -0.1260     *
## sp06_intercept  0.342000 0.0693  0.1990  0.4770     *
## sp06_env       -0.463000 0.0339 -0.4990 -0.3750     *
## sp06_env2      -0.329000 0.0655 -0.4550 -0.2050     *
## sp07_intercept  0.469000 0.0285  0.3950  0.4990     *
## sp07_env       -0.430000 0.0756 -0.4990 -0.2470     *
## sp07_env2      -0.441000 0.0451 -0.4980 -0.3320     *
## sp08_intercept  0.483000 0.0171  0.4370  0.5000     *
## sp08_env       -0.356000 0.0953 -0.4940 -0.1460     *
## sp08_env2      -0.460000 0.0338 -0.4990 -0.3780     *
## sp09_intercept  0.484000 0.0158  0.4390  0.5000     *
## sp09_env       -0.362000 0.0874 -0.4920 -0.1740     *
## sp09_env2      -0.478000 0.0211 -0.4990 -0.4220     *
## sp10_intercept  0.488000 0.0117  0.4550  0.5000     *
## sp10_env       -0.091000 0.0934 -0.2500  0.1240      
## sp10_env2      -0.484000 0.0144 -0.4990 -0.4470     *
## sp11_intercept  0.488000 0.0115  0.4590  0.5000     *
## sp11_env        0.198000 0.0886  0.0421  0.3900     *
## sp11_env2      -0.475000 0.0230 -0.4990 -0.4160     *
## sp12_intercept  0.487000 0.0130  0.4490  0.5000     *
## sp12_env        0.462000 0.0346  0.3710  0.4990     *
## sp12_env2      -0.482000 0.0176 -0.5000 -0.4360     *
## sp13_intercept  0.478000 0.0211  0.4210  0.4990     *
## sp13_env        0.423000 0.0643  0.2730  0.4980     *
## sp13_env2      -0.478000 0.0205 -0.4990 -0.4240     *
## sp14_intercept  0.475000 0.0238  0.4090  0.4990     *
## sp14_env        0.457000 0.0511  0.3140  0.4990     *
## sp14_env2      -0.454000 0.0381 -0.4980 -0.3600     *
## sp15_intercept  0.392000 0.0633  0.2600  0.4920     *
## sp15_env        0.466000 0.0318  0.3810  0.4990     *
## sp15_env2      -0.347000 0.0678 -0.4710 -0.2050     *
## sp16_intercept  0.276000 0.0727  0.1390  0.4220     *
## sp16_env        0.476000 0.0288  0.4080  0.5000     *
## sp16_env2      -0.259000 0.0701 -0.3960 -0.1260     *
## sp17_intercept  0.177000 0.0707  0.0404  0.3200     *
## sp17_env        0.478000 0.0211  0.4210  0.4990     *
## sp17_env2      -0.171000 0.0691 -0.3090 -0.0392     *
## sp18_intercept -0.000261 0.0741 -0.1480  0.1430      
## sp18_env        0.481000 0.0197  0.4280  0.4990     *
## sp18_env2       0.009230 0.0702 -0.1350  0.1450      
## sp19_intercept -0.092200 0.0739 -0.2300  0.0535      
## sp19_env        0.480000 0.0201  0.4240  0.4990     *
## sp19_env2       0.101000 0.0678 -0.0311  0.2390      
## sp20_intercept -0.277000 0.0738 -0.4190 -0.1340     *
## sp20_env        0.471000 0.0283  0.3960  0.4990     *
## sp20_env2       0.246000 0.0695  0.1060  0.3820     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.414, and the DIC is 77410.  Computation involved 2500 Gibbs steps, with a burnin of 500.  Dimension reduction was implemented with N = 11 and r = 8.
  Tau <- solve(mod_gjam1$parameters$sigMu)
  Tau_n = to_prec(Tau)
  
  #postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
  #postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
  
  
  #pH<-convert_to_m(postH)
  #pL<-convert_to_m(postL)
  
  #R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
  
  par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
  corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("Correlation cor(Y)")
  corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("R")
  corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("E matrix")
  # corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  # title("Tau")
  # corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  # title("R signif")
  corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
  title("True interactions")

data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod)

Environmental filtering + Facilitation dense 5 species

JSDM

mf5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mf5)
## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
#mf5$Rhat
cat(sprintf("Maximum Rhat value for Beta: %s\n", round(max(mf5$R$B),4)))
## Maximum Rhat value for Beta: 1.3081
cat(sprintf("Maximum Rhat value for Rho: %s\n", round(max(mf5$R$Rho),4)))
## Maximum Rhat value for Rho: NA
cat(sprintf("Maximum Rhat value for EnvRho: %s\n", round(max(mf5$R$EnvRho),4)))
## Maximum Rhat value for EnvRho: 1.1421
cat(sprintf("Maximum Rhat value for Tau: %s\n", round(max(mf5$R$Tau),4)))
## Warning in max(mf5$R$Tau): no non-missing arguments to max; returning -Inf
## Maximum Rhat value for Tau: -Inf
mf5$n.eff
## $B
##      [,1] [,2] [,3]
## [1,]   61   47   57
## [2,]   75  496   90
## [3,]  117  305  106
## [4,]   15   15   16
## [5,]  137  161  170
## 
## $Rho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1   44   44  678   61
## [2,]   44    1   88  199  334
## [3,]   44   88    1  282  110
## [4,]  678  199  282    1  171
## [5,]   61  334  110  171    1
## 
## $EnvRho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1  102   59   41   63
## [2,]  102    1  156   33  102
## [3,]   59  156    1   35  225
## [4,]   41   33   35    1   34
## [5,]   63  102  225   34    1
mf5$mcmc.info[1:7]
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000
data<-sim_data$FacDenseSp5
data <- list(
  Y = subset(data, select = -env),
  X = cbind(1, scale(poly(data$env, 2))),
  covx = cov(cbind(1, scale(poly(data$env, 2)))),
  K = 3,
  J = ncol(data) - 1,
  n = nrow(data),
  I = diag(ncol(data) - 1),
  df = ncol(data)
)

##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))

plot_cor_jsdm(mf5,data$Y,fac_inter[[4]])

Gjam

data<-sim_data$FacDenseSp5

#fit_gjam(data,5000,500,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
load_gjam(data,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])
## 
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env     113.0 61.5   13.9    237
## env2     67.7 31.3   19.6    134
## 
## Coefficient matrix B:
##             sp01   sp02   sp03  sp04   sp05
## intercept -1.300  1.020  2.990 1.570 -1.840
## env       -1.540 -1.640  0.114 2.570  1.760
## env2       0.558 -0.451 -1.120 0.308  0.299
## RMSPE      0.308  0.315  0.201 0.331  0.308
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025 CI_975 sig95
## sp01_intercept   -1.300 0.0975 -1.4900 -1.110     *
## sp01_env         -1.540 0.1450 -1.8300 -1.300     *
## sp01_env2         0.558 0.0909  0.3960  0.749     *
## sp02_intercept    1.020 0.1180  0.8440  1.300     *
## sp02_env         -1.640 0.0966 -1.8200 -1.450     *
## sp02_env2        -0.451 0.0802 -0.6250 -0.298     *
## sp03_intercept    2.990 0.1230  2.7500  3.230     *
## sp03_env          0.114 0.0787 -0.0114  0.287      
## sp03_env2        -1.120 0.0715 -1.2600 -0.977     *
## sp04_intercept    1.570 0.0908  1.4100  1.770     *
## sp04_env          2.570 0.1370  2.3300  2.870     *
## sp04_env2         0.308 0.0514  0.2070  0.410     *
## sp05_intercept   -1.840 0.1750 -2.1700 -1.510     *
## sp05_env          1.760 0.1520  1.4600  2.060     *
## sp05_env2         0.299 0.0823  0.1450  0.448     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 5 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.296, and the DIC is 80408.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gjfd5<-load_object("./gjam_models/gjam5f.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[4]])

Environmental filtering + Facilitation dense 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808
## Warning in max(mfd10$R$Tau): no non-missing arguments to max; returning -
## Inf
## Maximum Rhat value for Tau: -Inf
## $B
##       [,1] [,2] [,3]
##  [1,]   25   94   22
##  [2,]   27   26   25
##  [3,]  223  198  394
##  [4,]   24   27   23
##  [5,]  661 1000  896
##  [6,]   22   20   29
##  [7,]  275  464   65
##  [8,]  445  489  889
##  [9,]   16   18   19
## [10,]   19   33   16
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1   96  162   47  215   61  204  110  270    93
##  [2,]   96    1  200   34  843   55  883  202  223   443
##  [3,]  162  200    1  272 1000  801   53 1000  132   757
##  [4,]   47   34  272    1  387  126  281  247   83    72
##  [5,]  215  843 1000  387    1  233  210  467  590   357
##  [6,]   61   55  801  126  233    1   41   33   35    71
##  [7,]  204  883   53  281  210   41    1  184  531    55
##  [8,]  110  202 1000  247  467   33  184    1 1000   864
##  [9,]  270  223  132   83  590   35  531 1000    1   161
## [10,]   93  443  757   72  357   71   55  864  161     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1  109   29   42   43   25   42   29   15    15
##  [2,]  109    1   62   60   54   33   98   44   12    12
##  [3,]   29   62    1   26  580   27  200 1000   17    17
##  [4,]   42   60   26    1   30   34   37   32   24    22
##  [5,]   43   54  580   30    1   26  163  740   27    34
##  [6,]   25   33   27   34   26    1   46   88   26    22
##  [7,]   42   98  200   37  163   46    1  499   19    20
##  [8,]   29   44 1000   32  740   88  499    1   22    27
##  [9,]   15   12   17   24   27   26   19   22    1   105
## [10,]   15   12   17   22   34   22   20   27  105     1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

data<-sim_data$FacDenseSp10

#fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
load_gjam(data,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])
## 
## Sensitivity by predictor variables f:
##      Estimate    SE CI_025 CI_975
## env     160.0 35.60  101.0    240
## env2     18.5  4.88   10.4     29
## 
## Coefficient matrix B:
##             sp01   sp02   sp03   sp04    sp05   sp06   sp07   sp08  sp09
## intercept  0.212 -1.490 -0.877  2.410 -0.0951  2.400  0.497 -0.401 1.260
## env       -2.760 -2.660 -1.810 -3.010 -0.5910  1.160  1.330  1.030 3.990
## env2       1.070 -0.969 -0.644  0.319 -0.7020 -0.454 -0.256 -0.347 0.095
## RMSPE      0.274  0.375  0.406  0.305  0.4930  0.259  0.390  0.471 0.248
##             sp10
## intercept -0.893
## env        2.250
## env2      -0.411
## RMSPE      0.354
## 
## Coefficient matrix B:
##                Estimate     SE  CI_025  CI_975 sig95
## sp01_intercept   0.2120 0.0687  0.0833  0.3520     *
## sp01_env        -2.7600 0.1740 -3.1200 -2.4700     *
## sp01_env2        1.0700 0.1350  0.8190  1.3200     *
## sp02_intercept  -1.4900 0.2000 -1.9200 -1.1600     *
## sp02_env        -2.6600 0.2110 -3.0900 -2.3000     *
## sp02_env2       -0.9690 0.1130 -1.2400 -0.7780     *
## sp03_intercept  -0.8770 0.1070 -1.0900 -0.6900     *
## sp03_env        -1.8100 0.1470 -2.1200 -1.5700     *
## sp03_env2       -0.6440 0.1190 -0.9250 -0.4570     *
## sp04_intercept   2.4100 0.2040  2.0500  2.8000     *
## sp04_env        -3.0100 0.2470 -3.4600 -2.5900     *
## sp04_env2        0.3190 0.0696  0.1860  0.4540     *
## sp05_intercept  -0.0951 0.0554 -0.2050  0.0147      
## sp05_env        -0.5910 0.0628 -0.7140 -0.4700     *
## sp05_env2       -0.7020 0.0743 -0.8490 -0.5590     *
## sp06_intercept   2.4000 0.1260  2.1600  2.6400     *
## sp06_env         1.1600 0.0761  1.0100  1.3100     *
## sp06_env2       -0.4540 0.0649 -0.5840 -0.3310     *
## sp07_intercept   0.4970 0.0563  0.3840  0.6040     *
## sp07_env         1.3300 0.0797  1.1800  1.4900     *
## sp07_env2       -0.2560 0.0554 -0.3670 -0.1470     *
## sp08_intercept  -0.4010 0.0539 -0.5070 -0.2950     *
## sp08_env         1.0300 0.0877  0.8540  1.1900     *
## sp08_env2       -0.3470 0.0592 -0.4640 -0.2320     *
## sp09_intercept   1.2600 0.0922  1.0900  1.4500     *
## sp09_env         3.9900 0.3810  3.3800  4.6800     *
## sp09_env2        0.0950 0.0613 -0.0361  0.2050      
## sp10_intercept  -0.8930 0.0796 -1.0700 -0.7550     *
## sp10_env         2.2500 0.1280  2.0400  2.5200     *
## sp10_env2       -0.4110 0.0576 -0.5230 -0.2980     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 10 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.367, and the DIC is 157914.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#gjfd5<-load_object("./gjam_models/gjam10fd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

Environmental filtering + Facilitation dense 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508
## Warning in max(mfd20$R$Tau): no non-missing arguments to max; returning -
## Inf
## Maximum Rhat value for Tau: -Inf
## $B
##       [,1] [,2] [,3]
##  [1,]   13   57   12
##  [2,]   46   42   42
##  [3,]  301  176   90
##  [4,]   10   29   10
##  [5,]   46   37   35
##  [6,]   34   32   37
##  [7,]  536  181  102
##  [8,]  868 1000  388
##  [9,]  121  172  238
## [10,]  100  625  303
## [11,]  499 1000  121
## [12,] 1000  184  395
## [13,]  269  166  153
## [14,]   54  149  179
## [15,]  112  231  138
## [16,]  108  227   63
## [17,]   42  117   37
## [18,]   44   43   71
## [19,]   16   29   14
## [20,]   23   31   27
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1  353  436   46  381   42  169   83   93    24   171   130   160
##  [2,]  353    1  392  339   33   80  150  134   41    26   184   175   106
##  [3,]  436  392    1   51   59  101  171   59  297   569   166   359   111
##  [4,]   46  339   51    1 1000   90   35   68   27   118    78    50    40
##  [5,]  381   33   59 1000    1  269  156  867  283    48   130    61    46
##  [6,]   42   80  101   90  269    1   33   79  162   336    44    51   220
##  [7,]  169  150  171   35  156   33    1  303  143   127    97   190    66
##  [8,]   83  134   59   68  867   79  303    1  147   109   472   344    52
##  [9,]   93   41  297   27  283  162  143  147    1   292   314   137    62
## [10,]   24   26  569  118   48  336  127  109  292     1   367   362   180
## [11,]  171  184  166   78  130   44   97  472  314   367     1   326   229
## [12,]  130  175  359   50   61   51  190  344  137   362   326     1  1000
## [13,]  160  106  111   40   46  220   66   52   62   180   229  1000     1
## [14,]  243  134  169   37  160  166   73   44   36   131   100   182   173
## [15,]  113   82  622   46   37   51  136  145   72    46   273   159    85
## [16,]  555   52  111   58   47   64  285   64   22    28   151   525    36
## [17,]  360   84   62   79   73  569  191  106   45    52   120   102    29
## [18,]  129  128   86  191  621  197  226  202  257    53    32  1000   817
## [19,]  119   73  102  119  348   70  168   41   35   163    41    52    29
## [20,]  190  194  199  117  127  376   84  122  674    46   124    70    40
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]   243   113   555   360   129   119   190
##  [2,]   134    82    52    84   128    73   194
##  [3,]   169   622   111    62    86   102   199
##  [4,]    37    46    58    79   191   119   117
##  [5,]   160    37    47    73   621   348   127
##  [6,]   166    51    64   569   197    70   376
##  [7,]    73   136   285   191   226   168    84
##  [8,]    44   145    64   106   202    41   122
##  [9,]    36    72    22    45   257    35   674
## [10,]   131    46    28    52    53   163    46
## [11,]   100   273   151   120    32    41   124
## [12,]   182   159   525   102  1000    52    70
## [13,]   173    85    36    29   817    29    40
## [14,]     1    90    96   127   116   154    43
## [15,]    90     1   127   137    82    28   237
## [16,]    96   127     1   279   122    95   192
## [17,]   127   137   279     1   136   151   227
## [18,]   116    82   122   136     1    90   236
## [19,]   154    28    95   151    90     1    57
## [20,]    43   237   192   227   236    57     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1   16   18   17   18   15   14   15   15    17    18    17    15
##  [2,]   16    1  105   12   62   36   45   71   75   190   103   139    74
##  [3,]   18  105    1   17   87   57  112  130  220    93   266   109   139
##  [4,]   17   12   17    1   11   17   12   12   12    12    15    13    12
##  [5,]   18   62   87   11    1   43   62   92  135   162   193   102   124
##  [6,]   15   36   57   17   43    1   48   39   39    32    47    33    35
##  [7,]   14   45  112   12   62   48    1  551  219   209   874   259   208
##  [8,]   15   71  130   12   92   39  551    1  755   425  1000   409   416
##  [9,]   15   75  220   12  135   39  219  755    1   407  1000   501  1000
## [10,]   17  190   93   12  162   32  209  425  407     1   668  1000   340
## [11,]   18  103  266   15  193   47  874 1000 1000   668     1   708  1000
## [12,]   17  139  109   13  102   33  259  409  501  1000   708     1   554
## [13,]   15   74  139   12  124   35  208  416 1000   340  1000   554     1
## [14,]   15   78   76   12  133   36  314  516  264   557   591   322   306
## [15,]   14   72  147   12   83   32  133  231 1000   222   611   434  1000
## [16,]   15   26  138   15   40   93   84   70   78    53   119    59    78
## [17,]   80   89   71   23   76   41   35   42   64    55    57    57    50
## [18,]   19   44   92   16   39   44  230  151  136   127   223   241   132
## [19,]   32   22   17   38   21   23   17   17   18    22    22    23    18
## [20,]   21   78   48   17   40   21   27   31   34    39    38    34    30
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]    15    14    15    80    19    32    21
##  [2,]    78    72    26    89    44    22    78
##  [3,]    76   147   138    71    92    17    48
##  [4,]    12    12    15    23    16    38    17
##  [5,]   133    83    40    76    39    21    40
##  [6,]    36    32    93    41    44    23    21
##  [7,]   314   133    84    35   230    17    27
##  [8,]   516   231    70    42   151    17    31
##  [9,]   264  1000    78    64   136    18    34
## [10,]   557   222    53    55   127    22    39
## [11,]   591   611   119    57   223    22    38
## [12,]   322   434    59    57   241    23    34
## [13,]   306  1000    78    50   132    18    30
## [14,]     1   248    61    40   106    21    29
## [15,]   248     1    63    54   124    18    30
## [16,]    61    63     1    29   124    14    24
## [17,]    40    54    29     1    49    59   112
## [18,]   106   124   124    49     1    23    25
## [19,]    21    18    14    59    23     1    42
## [20,]    29    30    24   112    25    42     1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Gjam

data<-sim_data$FacDenseSp20

fit_gjam(data,5000,500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
## 
## Observations and responses:
## [1] 500  20
## ===========================================================================
## expanding covariance chains
## ===========================================================================
## Sensitivity by predictor variables f:
##      Estimate   SE CI_025 CI_975
## env       307 54.9    217    427
## env2      187 34.7    122    262
## 
## Coefficient matrix B:
##             sp01     sp02   sp03   sp04    sp05   sp06    sp07   sp08
## intercept -0.480 -0.99800 -0.423 -0.180 -0.0128  0.385  0.0559  0.244
## env       -1.980 -1.74000 -1.860 -2.680 -1.3400 -1.960 -0.9680 -0.722
## env2       0.308  0.00803 -0.255 -0.555 -0.4870 -0.893 -0.6740 -0.899
## RMSPE      0.321  0.35800  0.370  0.296  0.4070  0.291  0.4510  0.445
##             sp09   sp10   sp11   sp12   sp13   sp14    sp15   sp16    sp17
## intercept  0.640  0.712  0.317  0.262  0.579  0.250  0.0141 -0.209 -0.0225
## env       -0.914 -0.402  0.105  0.469  0.932  1.050  1.5200  2.100  2.1400
## env2      -1.250 -1.120 -0.865 -0.937 -0.874 -0.957 -1.0500 -0.612 -0.1930
## RMSPE      0.383  0.413  0.471  0.453  0.395  0.407  0.3710  0.331  0.3260
##             sp18  sp19   sp20
## intercept 0.0822 0.125 -0.385
## env       2.4300 2.800  2.080
## env2      0.2020 0.794  0.260
## RMSPE     0.2980 0.278  0.317
## 
## Coefficient matrix B:
##                Estimate     SE   CI_025  CI_975 sig95
## sp01_intercept -0.48000 0.0709 -0.61400 -0.3370     *
## sp01_env       -1.98000 0.1840 -2.28000 -1.6500     *
## sp01_env2       0.30800 0.0638  0.18000  0.4320     *
## sp02_intercept -0.99800 0.1380 -1.27000 -0.7540     *
## sp02_env       -1.74000 0.1240 -1.99000 -1.5200     *
## sp02_env2       0.00803 0.0846 -0.16100  0.1510      
## sp03_intercept -0.42300 0.1680 -0.75800 -0.2050     *
## sp03_env       -1.86000 0.1500 -2.12000 -1.5500     *
## sp03_env2      -0.25500 0.1150 -0.45600 -0.0566     *
## sp04_intercept -0.18000 0.0555 -0.28400 -0.0712     *
## sp04_env       -2.68000 0.1380 -2.95000 -2.4200     *
## sp04_env2      -0.55500 0.0641 -0.67700 -0.4330     *
## sp05_intercept -0.01280 0.0903 -0.17800  0.1530      
## sp05_env       -1.34000 0.0753 -1.49000 -1.2000     *
## sp05_env2      -0.48700 0.1010 -0.65700 -0.2960     *
## sp06_intercept  0.38500 0.0585  0.27000  0.4990     *
## sp06_env       -1.96000 0.1200 -2.21000 -1.7400     *
## sp06_env2      -0.89300 0.0928 -1.06000 -0.6990     *
## sp07_intercept  0.05590 0.0662 -0.07310  0.1860      
## sp07_env       -0.96800 0.0724 -1.12000 -0.8310     *
## sp07_env2      -0.67400 0.0632 -0.80400 -0.5520     *
## sp08_intercept  0.24400 0.0603  0.13200  0.3630     *
## sp08_env       -0.72200 0.0649 -0.84700 -0.5930     *
## sp08_env2      -0.89900 0.0723 -1.05000 -0.7690     *
## sp09_intercept  0.64000 0.0555  0.53100  0.7490     *
## sp09_env       -0.91400 0.0596 -1.03000 -0.7970     *
## sp09_env2      -1.25000 0.0745 -1.40000 -1.1100     *
## sp10_intercept  0.71200 0.0581  0.59800  0.8230     *
## sp10_env       -0.40200 0.0555 -0.51100 -0.2930     *
## sp10_env2      -1.12000 0.0838 -1.27000 -0.9570     *
## sp11_intercept  0.31700 0.0807  0.13400  0.4750     *
## sp11_env        0.10500 0.0607 -0.01620  0.2190      
## sp11_env2      -0.86500 0.0678 -0.99500 -0.7330     *
## sp12_intercept  0.26200 0.0541  0.15600  0.3680     *
## sp12_env        0.46900 0.0547  0.36200  0.5740     *
## sp12_env2      -0.93700 0.0660 -1.06000 -0.8050     *
## sp13_intercept  0.57900 0.0808  0.41100  0.7210     *
## sp13_env        0.93200 0.0628  0.81000  1.0500     *
## sp13_env2      -0.87400 0.0659 -1.01000 -0.7510     *
## sp14_intercept  0.25000 0.0625  0.13100  0.3720     *
## sp14_env        1.05000 0.0695  0.90700  1.1800     *
## sp14_env2      -0.95700 0.1110 -1.14000 -0.7350     *
## sp15_intercept  0.01410 0.0589 -0.09820  0.1310      
## sp15_env        1.52000 0.0761  1.38000  1.6700     *
## sp15_env2      -1.05000 0.0809 -1.21000 -0.9040     *
## sp16_intercept -0.20900 0.0883 -0.37100 -0.0447     *
## sp16_env        2.10000 0.1020  1.91000  2.3100     *
## sp16_env2      -0.61200 0.0559 -0.72300 -0.5050     *
## sp17_intercept -0.02250 0.0508 -0.12200  0.0783      
## sp17_env        2.14000 0.1300  1.92000  2.4000     *
## sp17_env2      -0.19300 0.0816 -0.34100 -0.0369     *
## sp18_intercept  0.08220 0.0547 -0.02540  0.1870      
## sp18_env        2.43000 0.1090  2.22000  2.6500     *
## sp18_env2       0.20200 0.0628  0.07750  0.3160     *
## sp19_intercept  0.12500 0.0624  0.00441  0.2470     *
## sp19_env        2.80000 0.1190  2.57000  3.0300     *
## sp19_env2       0.79400 0.0876  0.63200  0.9540     *
## sp20_intercept -0.38500 0.0727 -0.51200 -0.2370     *
## sp20_env        2.08000 0.0929  1.90000  2.2600     *
## sp20_env2       0.26000 0.0780  0.10300  0.3960     *
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
## Coefficient matrix B, standardized for X and W:
## NULL
## 
## Last column indicates if 95% posterior distribution contains zero.
## 
##  Design Table
##      env env2
## VIF    1    1
## env2   0   NA
## 
##  Sample contains n = 500 observations on S = 20 response variables, and 2 predictors.  Data types (typeNames) include PA. There are 0 missing values in X and 0 missing values in Y. The RMSPE is 0.374, and the DIC is 395024.  Computation involved 5000 Gibbs steps, with a burnin of 500.

#load_gjam(data,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")

#to check posterior density of s in Sigma 
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))

HMSC

data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"L_d",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)

hm_inter(hm_mod, interact = fac_inter[[6]])

Environmental filtering + Facilitation sparse 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048
## Warning in max(mfs5$R$Tau): no non-missing arguments to max; returning -Inf
## Maximum Rhat value for Tau: -Inf
## $B
##      [,1] [,2] [,3]
## [1,]   17   44   17
## [2,]  169  469  127
## [3,] 1000 1000  500
## [4,]   20   25   20
## [5,]   13   66   14
## 
## $Rho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1  111   92  123  819
## [2,]  111    1  283   40  143
## [3,]   92  283    1  205   28
## [4,]  123   40  205    1  102
## [5,]  819  143   28  102    1
## 
## $EnvRho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1   20   19   29   13
## [2,]   20    1  252   20   18
## [3,]   19  252    1   21   17
## [4,]   29   20   21    1   32
## [5,]   13   18   17   32    1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Environmental filtering + Facilitation sparse 10 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746
## Warning in max(mfs10$R$Tau): no non-missing arguments to max; returning -
## Inf
## Maximum Rhat value for Tau: -Inf
## $B
##       [,1] [,2] [,3]
##  [1,]   52   49   55
##  [2,]   11   17   10
##  [3,]   19   20   24
##  [4,] 1000  374  266
##  [5,] 1000  567  333
##  [6,]  262  286  296
##  [7,]  839  602 1000
##  [8,]  129  212   94
##  [9,]  462  346  407
## [10,]   49  251   47
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1  122   94  384  296  316  114  140  292   337
##  [2,]  122    1  558  143  481  155   81  262  144   148
##  [3,]   94  558    1   89  355  800 1000  317  214   104
##  [4,]  384  143   89    1  172  290  191  213  543   438
##  [5,]  296  481  355  172    1  273  110  421  192   116
##  [6,]  316  155  800  290  273    1  621  178  401   101
##  [7,]  114   81 1000  191  110  621    1  179  140    53
##  [8,]  140  262  317  213  421  178  179    1  630   181
##  [9,]  292  144  214  543  192  401  140  630    1   213
## [10,]  337  148  104  438  116  101   53  181  213     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1   15   31  147  111  128  138   80  119    41
##  [2,]   15    1   37   13   14   15   12   22   11    15
##  [3,]   31   37    1   26   25   30   25   30   23    97
##  [4,]  147   13   26    1  449  363 1000   78 1000    48
##  [5,]  111   14   25  449    1  615  361   97  630    70
##  [6,]  128   15   30  363  615    1  696  109  358    50
##  [7,]  138   12   25 1000  361  696    1  110  488    48
##  [8,]   80   22   30   78   97  109  110    1  153    34
##  [9,]  119   11   23 1000  630  358  488  153    1    68
## [10,]   41   15   97   48   70   50   48   34   68     1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Environmental filtering + Facilitation sparse 20 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826
## Warning in max(mfs20$R$Tau): no non-missing arguments to max; returning -
## Inf
## Maximum Rhat value for Tau: -Inf
## $B
##       [,1] [,2] [,3]
##  [1,]   24   21   32
##  [2,]   22   59   19
##  [3,]   27  108   32
##  [4,]   22  251   24
##  [5,]  227   23  216
##  [6,]  273  395  140
##  [7,]  351 1000  250
##  [8,]  177  104  232
##  [9,]  206  144   69
## [10,]  384  290  280
## [11,]  128  404 1000
## [12,] 1000  570  505
## [13,]  887  747  270
## [14,]  170 1000 1000
## [15,]  175   51  217
## [16,]   43  147   29
## [17,]   10   29   13
## [18,]   22   26   16
## [19,]   12   84   13
## [20,]   23   33   25
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1  422   57   96   17   43   24   74   81    32   502   101    59
##  [2,]  422    1   91  212   54  128  110  167   93    75    89   128    87
##  [3,]   57   91    1  172   94  434  162  154  111   298   118    29    33
##  [4,]   96  212  172    1   92  399  138   52  123    41    74    65    43
##  [5,]   17   54   94   92    1  164   95  612  312    19    71    18    95
##  [6,]   43  128  434  399  164    1  152  162  180   349   118    37    80
##  [7,]   24  110  162  138   95  152    1  379  185   129   435    48   101
##  [8,]   74  167  154   52  612  162  379    1  105   519   122   477   204
##  [9,]   81   93  111  123  312  180  185  105    1   112    61   263   190
## [10,]   32   75  298   41   19  349  129  519  112     1   247   184   102
## [11,]  502   89  118   74   71  118  435  122   61   247     1   120   301
## [12,]  101  128   29   65   18   37   48  477  263   184   120     1    92
## [13,]   59   87   33   43   95   80  101  204  190   102   301    92     1
## [14,]  150   38   72  603   57  373   83  739   91   354    80   100   193
## [15,]   45   80   61  273   90   38   24   82   39   153    99    95   137
## [16,]   92  111   36   23   43   83  107   95  198    34    55   302    66
## [17,]   97   28   53  253  137  244  115   38  170    26    25   104   118
## [18,]  107   38  109  107  103   72   57   73   98   121    44    47    27
## [19,]   46  150  110  517   52  216   66   81   42    23   137    80    54
## [20,]   51   66  339  154   64   75   85   98   68    98    82    40    31
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]   150    45    92    97   107    46    51
##  [2,]    38    80   111    28    38   150    66
##  [3,]    72    61    36    53   109   110   339
##  [4,]   603   273    23   253   107   517   154
##  [5,]    57    90    43   137   103    52    64
##  [6,]   373    38    83   244    72   216    75
##  [7,]    83    24   107   115    57    66    85
##  [8,]   739    82    95    38    73    81    98
##  [9,]    91    39   198   170    98    42    68
## [10,]   354   153    34    26   121    23    98
## [11,]    80    99    55    25    44   137    82
## [12,]   100    95   302   104    47    80    40
## [13,]   193   137    66   118    27    54    31
## [14,]     1    69   337    73    53    44    61
## [15,]    69     1    88    19    33    36    34
## [16,]   337    88     1   210    83    48    51
## [17,]    73    19   210     1    51   125    51
## [18,]    53    33    83    51     1   123    37
## [19,]    44    36    48   125   123     1    37
## [20,]    61    34    51    51    37    37     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13]
##  [1,]    1   42   28   93   37   48   50   58   36    53    87    66    47
##  [2,]   42    1   15  128   56   25   37   24   36    35    29    28    29
##  [3,]   28   15    1   29   57   32   42   62   37    51    84    56    50
##  [4,]   93  128   29    1   25   32   26   32   37    50    39    31    33
##  [5,]   37   56   57   25    1  137  521  146  787   251   123   139   141
##  [6,]   48   25   32   32  137    1  536  243  205   335   240   303   169
##  [7,]   50   37   42   26  521  536    1  189  402   349   214   254   319
##  [8,]   58   24   62   32  146  243  189    1  251   258  1000  1000   283
##  [9,]   36   36   37   37  787  205  402  251    1   705   162   189   266
## [10,]   53   35   51   50  251  335  349  258  705     1   244   246  1000
## [11,]   87   29   84   39  123  240  214 1000  162   244     1  1000   367
## [12,]   66   28   56   31  139  303  254 1000  189   246  1000     1   324
## [13,]   47   29   50   33  141  169  319  283  266  1000   367   324     1
## [14,]   36   28   50   43  160  187  217  312  412  1000   276   246   749
## [15,]   55   31   48   36   89  301  132  550  140   294  1000  1000   305
## [16,]   72   22   40   57   26   43   37   50   46    76    73    56    58
## [17,]   17   47   24   37   15   15   14   20   15    17    26    23    15
## [18,]   76   74   21   51   23   30   21   33   22    29    51    36    23
## [19,]  108  155   40   28   32   17   19   19   24    22    21    18    18
## [20,]   41  182  143   38   68   42   39   31   55    57    34    33    38
##       [,14] [,15] [,16] [,17] [,18] [,19] [,20]
##  [1,]    36    55    72    17    76   108    41
##  [2,]    28    31    22    47    74   155   182
##  [3,]    50    48    40    24    21    40   143
##  [4,]    43    36    57    37    51    28    38
##  [5,]   160    89    26    15    23    32    68
##  [6,]   187   301    43    15    30    17    42
##  [7,]   217   132    37    14    21    19    39
##  [8,]   312   550    50    20    33    19    31
##  [9,]   412   140    46    15    22    24    55
## [10,]  1000   294    76    17    29    22    57
## [11,]   276  1000    73    26    51    21    34
## [12,]   246  1000    56    23    36    18    33
## [13,]   749   305    58    15    23    18    38
## [14,]     1   255    75    16    24    19    39
## [15,]   255     1   125    19    38    15    36
## [16,]    75   125     1    13    48    14    30
## [17,]    16    19    13     1   103    14    14
## [18,]    24    38    48   103     1    18    25
## [19,]    19    15    14    14    18     1    21
## [20,]    39    36    30    14    25    21     1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Environmental filtering + Competition dense 5 species

JSDM

## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## Successful convergence based on Rhat values (all < 1.1).
## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313
## Warning in max(mcmpd5$R$Tau): no non-missing arguments to max; returning -
## Inf
## Maximum Rhat value for Tau: -Inf
## $B
##      [,1] [,2] [,3]
## [1,] 1000 1000  389
## [2,]  210  355  355
## [3,]  320 1000  272
## [4,] 1000 1000 1000
## [5,] 1000  463   95
## 
## $Rho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1  252  209  114   59
## [2,]  252    1  132  483   52
## [3,]  209  132    1  611   96
## [4,]  114  483  611    1  162
## [5,]   59   52   96  162    1
## 
## $EnvRho
##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1  556  455  431  207
## [2,]  556    1  795  661  120
## [3,]  455  795    1 1000  130
## [4,]  431  661 1000    1  100
## [5,]  207  120  130  100    1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Environmental filtering + Competition dense 10 species

## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho 
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
## 
## For each of 5 chains:
## Adaptation:            250000 iterations (possibly insufficient)
## Burn-in:               0 iterations
## Thin rate:             100 iterations
## Total chain length:    270000 iterations
## Posterior sample size: 200 draws
## 
## **WARNING** Rhat values indicate convergence failure.
## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771
## Warning in max(mcmpd10$R$Tau): no non-missing arguments to max; returning -
## Inf
## Maximum Rhat value for Tau: -Inf
## $B
##       [,1] [,2] [,3]
##  [1,]   24   38   25
##  [2,]  802  157   93
##  [3,]  833  870 1000
##  [4,]  573  462  566
##  [5,]  994  601  648
##  [6,]  708 1000  329
##  [7,]  420  240  320
##  [8,]   29   28   33
##  [9,]  123  174  114
## [10,]  206  117  100
## 
## $Rho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1   49  197   83   87  125   89  176  184   124
##  [2,]   49    1   40   82  315   39   73  187   71   174
##  [3,]  197   40    1  387  377  273  177  209  857   269
##  [4,]   83   82  387    1   81  136 1000   50  131   442
##  [5,]   87  315  377   81    1 1000  118  169   59   219
##  [6,]  125   39  273  136 1000    1 1000  144  242   144
##  [7,]   89   73  177 1000  118 1000    1   55  216    41
##  [8,]  176  187  209   50  169  144   55    1  381   155
##  [9,]  184   71  857  131   59  242  216  381    1   175
## [10,]  124  174  269  442  219  144   41  155  175     1
## 
## $EnvRho
##       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
##  [1,]    1   39   33   59   59   34   24  134   31    59
##  [2,]   39    1   96   83   84  108  170  178   92   306
##  [3,]   33   96    1  509  591  612  105   37  105    96
##  [4,]   59   83  509    1 1000  326  105   44  114    97
##  [5,]   59   84  591 1000    1  385  119   39  181    99
##  [6,]   34  108  612  326  385    1  256   59  112   172
##  [7,]   24  170  105  105  119  256    1  113   67   194
##  [8,]  134  178   37   44   39   59  113    1   69   222
##  [9,]   31   92  105  114  181  112   67   69    1   123
## [10,]   59  306   96   97   99  172  194  222  123     1
## $n.chains
## [1] 5
## 
## $n.adapt
## [1] 250000 250000 250000 250000 250000
## 
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
## 
## $n.iter
## [1] 20000
## 
## $n.burnin
## [1] 0
## 
## $n.thin
## [1] 100
## 
## $n.samples
## [1] 1000

Environmental filtering + Competition dense 20 species

Environmental filtering + Competition sparse 5 species

Environmental filtering + Competition dense 10 species

Environmental filtering + Competition dense 5 species

Graph

Visualization for true interactions